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Abstract 

The  state  estimation  problem  for  bilinear  stochastic  systems 
evolving  on  compact  Lie  groups  and  homogeneous  spaces  is 
considered . The  problem  is  motivated  by  some  applications 
involving  rotational  processes  in  three  dimensions. 

The  theory  of  harmonic  analysis  on  compact  Lie  groups  is 
used  to  define  assumed  density  approximations  which  re- 
sult in  implementable  suboptimal  estimators  for  the  state  of 
the  bilinear  system.  The  results  of  Monte  Carlo  simulations 
are  reported;  these  indicate  that  simple  filters  designed  by 
these  techniques  perform  well  as  compared  to  other  filters. 
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Fourier  series  analysis  has  been  applied  in  several  recent  studies 
[l]  - [4]  to  estimation  problems  for  stochastic  processes  evolving  on  the 
circle  S^.  Willsky  [4]  used  Fourier  series  methods  to  define  "assumed 
density"  approximations  for  certain  phase  tracking  and  demodulation 
problems.  In  fact,  a system  designed  using  these  techniques  performed 
better  than  other  estimators.  Including  an  optimal  phase-lock  loop. 

In  this  paper  we  study  bilinear  systems  evolving  on  compact  He 
groups  or  homogeneous  spaces  [29].  The  optimal  estimator  is  . in 
general  infinite  dimensional  [7] , and  our  approach  to  the  design  of 
suboptimal  estimators  is  a generalization  of  that  of  Willsky  [4],  whose  work  is 
reviewed  briefly  in  Section  III.  The  basic  approach  involves  the  definition 
of  an  "assumed  density"  form  for  the  conditional  density  of  the  system 
i state  at  time  t given  observations  up  to  time  t.  These  densities  are  defined 
via  the  techniques  of  harmonic  analysis  on  compact  Lie  groups  [5], [6] 

(which  generalize  the  Fourier  series  on  the  Lie  group  S^) . Our  method  differs 
frpm  most  previous  assumed  density  approximations  in  that  our  approximation 
is  defined  on  the  appropriate  compact  manifold  (as  opposed  to  the  usual 
Gaussian  approximations,  for  example,  which  are  defined  on  IR^  [7]). 

This  method  also  avoids  the  problem  of  merely  truncating  higher  order  terms 
in  a harmonic  expansion;  as  pointed  out  by  Lo  [18]  and  Willsky  [4],  such 
higher  order  terms  will  not  be  negligible,  especially  if  the  filter  is 
performing  well.  For  an  alternative  approach  to  discrete-time  estimation 
problems  on  Lie  groups  and  homogeneous  spaces,  see  the  work  of  Lo  and 
Eshleman  [18-20],  who  use  exponential  Fourier  densities  to  avoid  the  trun- 
cation problem. 

In  Section  II  we  review  some  general  properties  of  stochastic 


bilinear  systems  and  discuss  the  estimation  problem  for  systems  evolving 
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on  compact  Lie  groups  and  homogeneous  spaces*  Section  111  contains  the 

application  of  the  technique  to  systems  evolving  on  s”,  while  Section  IV 

contains  the  application  to  systems  on  SOCn).  Results  of  Monte  Carlo 

2 

simulations  of  the  S estimator  are  presented  in  Section  Finally,  the 
Appendix  reviews  some  necessary  concepts  from  the  theory  of  harmonic  analy- 
sis on  compact  Lie  groups.  Some  of  the  concepts  of  this  paper  were  intro- 
duced in  [21]  and  [22],  hut  no  simulation  results  were  presented. 
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II,  ESTIMATION  FOR  STOCHASTIC  BILINEAR  SYSTEMS 

The  basic  stochastic  bilinear  system  Cor  linear  system  with  state- 
dependent  noise)  considered  here  is  described  by  the  Ito  stochastic  diffe- 
rential equation  [4],  [9-17],  [21],  [25] 


N 


N 


dxCt)  “{[A^  + Ja  z Q,.(t)A  A ]dt  + Z A.dw  (t)}xCt)  (1) 
l,j=l  -*  i”l 

where  x is  an  n-vector  or  an  n x n matrix,,  the  A^  are  n x n matrices, 

is  the  Ci»j)th  element  of  Q,  and  v is  a Brownian  motion  (Wiener)  process 

mlnCt,s) 

with,  strength  QCt)  such  that  E[w(t)w’(s)]  “ / Q(T)dT.  Following 

0 

the  notation  of  [8-11],  we  define  ,2’=  {Aq,Aj^,  ..  to  be  the  smallest 

Lie  algebra  containing  these  matrices.  The  corresponding  connected  matrix 
Lie  group  is  denoted  by  G ■>  {exp^}-.  Then,  if  x is  an  n x n matrix  and 
xCtQ)EG,  the  solution  x(t)  of  (1)  evolves  on  G (l.e.,  x(t)eG  for  all  t > 0) 
in  the  mean-square  sense  and  almost  surely  [15-17].  If  x is  an  n-vector, 
then  the  solution  of  (1)  evolves  on  the  homogeneous  space  G'xCtg). 

Associated  with  the  Ito  equation  (1)  is  a sequence  of  equations  for 
the  powers  of  the  state  xCt)  (see  Brockett  [9],  [10]).  If  N(n,p)  denotes 


[pl 


the  binomial  coefficient  given  an  n-vector  x,  we  define  x 

to  be  the  N(n,p) -vector  with  components  equal  to  the  monomials  (homogeneous 
polynomials)  of  degree  p in  x^,.,.,x^,  the  components  of  x,  scaled  so  that  . 

11*11  ^ II  • Given  an  m x n matrix  A,  we  denote  by  A^^^  the  unique 

matrix  which  verifies 


, - Ax  ^ y'’’'  - a'p’  .'p' 


(2) 
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can  be  Interpreted  as  a linear  operator  on  symmetric  tensors  of  degree  p 
[9],  and  is  known  as  the  symmetrized  Kronecker  p^^  power  of  A [26],  It  is 
clear  that  if  x satisfies  the  linear  differential  equation 


x(t)  = Ax(t) 


then  also  satisfies  a linear  differential  equation 

x^^(t)  = A^p^x^P^(t)  . (4) 

We  regard  this  as  the  definition  of  Ar  which  is  the  infinitesimal  version  of 
[ ] LpJ 

A^*^  . In  fact,  Ar  , can  be  easily  computed  from  A [25]. 

LPJ  [p] 

It  can  easily  be  shown  that  if  x satisfies  (1),  then  satisfies 

the  Ito  equation 

dx^P^(t)={A-  +iEQ,Jt)A.  A 1x^P^(t)dt+  eA.  x^P^(t)dwJt)  (5) 

°Ep]  ^[p]  ^[p]  1=1  V]  ^ 

In  addition,  if  the  nxn  matrix  X satisfies  (1),  it  is  easy  to  show  that 
also  satisfies  (5) . As  we  shall  see  later  in  the  section,  this  sequence  of 
•equations  is  a valuable  tool  in  the  study  of  state  estimation. 

The  observation  model  considered  in  this  paper  consists  of  linear 
observations  of  the  state  corrupted  by  additive  white  noise,  or 

dz(t)  = L(x(t))  dt  dv(t)  (6) 

where  L is  a linear  operator  and  v is  a Wiener  process.  This  bilinear  system- 

linear  observation  model  is  useful  in  the  study  of  certain  practical  problems, 

2 

such  as  the  S satellite  tracking  and  S0(3)  rigid  body  orientation  estimation 
problems  discussed  in  [13,  Ch.  4]  and  [21,  Sec.  IV]. 

The  remainder  of  the  paper  is  devoted  to  the  study  of  the  estimation 
problem  for  two  classes  of  systems  of  the  form  (1) , (6) , which  are  motivated 
by  the  aforementioned  examples.  Tl\e  first  system  consists  of  the  bilinear 
state  equation  , 
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(7) 


N N 

dX(t)  = [A  4-  ^ r Q..(t)A  A ]X(t)clt+  EA  X(t)dw  (t) 

^ ij=l  J i=l  ^ 

with  linear  measurements 

dZj(t)  = X(t)h(t)dt+-  R^(t)dv(t)  (8) 

where  X(t)  and  {A^}  are  nxn  matrices,  Zj(t)  is  a p-vector,  w is  a Wiener 
process  with  strength  Q(t)  s 0,  v is  a standard  Wiener  process  independent 
of  w,  and  R > 0.  More  general  linear  measurerhents  can  obviously  be  con- 
sidered, but  for  simplicity  of  notation  we  restrict  our  attention  to  (8), 
which  arises  in  the  star  tracking  example  of  [13,  Ch.  4].  We  will  make 
the  crucial  assumption  that  the  Lie  group  G = {expi?}g  is  compact. 

The  second  system  consists  of  the  bilinear  state  equation  Cl)  with 
linear  measurements 

dz2Ct)  = HCt)x(t)dt  + R^Ct)dvCt)  (9) 

where  x(,t)  is  an  n-vector,  A^.  are  nxn  matrices,  and  Z2,  v and  w are  as 
above.  It  will  be  assumed  that  x evolves  on  a compact  homogeneous  space 
[8],  [13],  [29] — i.e.,  the  solution  of  Cl)  is 

xCt)  = X(t)x(0)  (10) 

where  X satisfies  C7)  with  XCO)  = I and  evolves  on  the  compact  Lie  group 
G ■ {exp5’}g. 

It  is  shown  in  [13]  that,  by  a linear  change  of  basis  on  the  state 
space,  C7)  and  Cl)  can  be  transformed  into  equations  which  evolve  on  the 
compact  special  orthogonal  group  SOCn)  *=  {XcIr" '^[X’X  = I,  det  X«+l}  and 
the  compact  homogeneous  space  S {xeIR  |x*x»'l}  Cthe  (n—1)— ‘Sphere) , rcspec— 

tlvely.  Hence,  we  need  only  consider  systems  evolving  on  SOCn)  and  S 
(this  is  equivalent  to  the  assumption  that  are  skew-symmetric 

1131), 
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ViU.iJJ 


The  estimation  criterion  which  will  be  used  for  these  two  problems 
is  the  constrained  least-squares  estimator,  which  is  analogous  to  the  cri- 
terion used  in  [1],  [4],  and  [21]  for  the  phase  estimation  problem.  That 
is,  for  C7)-C8)  we  wish  to  find  XCtjt)  which  minimizes  the  conditional  error 
covariance 

= E[tr{(XCt)  - XCt|t))’CX(t)  - X(tlt))}l2j]  Cll) 

subject  to  the  SOCn)  constraint  X(t| t) ’XCt j t)  = I,  where  the  notation  (11) 
denotes  the  conditional  expectation  given  the  a-field  o{z^}  generated  by 
■ e observed  process  z^  = (z^j^Cs),  0 < s < t}  up  to  time  t.  For  (!)»  (9) 
seek  xCt| t)  which  minimizes 

= E[CxCt)  - iCtlt))’  (xCt)  - i(tlt))Iz^]  (12) 

subject  to  the  S^~  constraint  HxCtlt)]]  = x(t|  t) 'x(tl  t)  = 1. 

It  is  easily  shown  [13]  that  the  optimal  estimates  are,  respectively 


X'ct  |t) =±x{t  1 t)[x{t  jt)  ’ X(t  1 1)]"^  , 

(13) 

_ x(t[tL. 
p(t|t)|| 

(14) 

the  conditional  expectation  is  denoted  by 

the  equivalent  notations 

xCtlt)  = E[x(.t)|22J  = E'^IxCt)]  , 

CIS) 

The  sign  in  Cl3)  is  chosen  to  Insure  that  det  X(tlt)  = +1  [23],  Thus  in 
both  cases  we  must  compute  the  conditional  expectation  of  the  state  (X(t)  or 
x(t))  given  the  past  obseirvatlons  (z^  = {zj^(s),  0 < s < t}  or  (z^Cs), 

0 < s < t}). 
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The  equations  for  computing  the  conditional  expectation  can  be  derived 
from  the  general  nonlinear  filtering  equation  [7]  and  the  moment  equation 
(5).  The  resultant  equations  for  the  SO(n)  system  (7)-(8)  jice 

dEWP^(tJ]=  [(A  E Q.,.(t)A  A )@l]E^[X^P^(t)]dt 

°[p]  ij=l  %]  J[p] 

+ [E‘[X^P^(t)h'{t)X(t)]-E‘[X^P^(t)]h’(t)E^[X(t)]  }R”^t)dVj(t)  (16) 

dv^(t)  = d2j(t)  - X(t|t)h(t)dt  (17) 

where  ^ denotes  Kronecker  product  and  X^^^  is  the  vector  containing  the 
elements  of  the  matrix  X^^"^  in  lexicographic  order  [26],  [32,  p.64].  For  the 
S**  ^ system  (1),  (9),  we  have 

dE‘[x^P^t)]  = [A-  +1  E Q.i(t)A.  A ]E‘[x*^P^(t)]dt 

%]  ij=l  Ip]  ^[p] 

+ [E^[x^P^  (t)x-  (t)]-E^[x^P^  (t)]E^[x’  (t)]  }H*  {t)R"^  (t)dv2  (t)  (18) 

dV2(t)  = dz^ft)  - H(t)x(t|t)dt  . (19) 

The  structure  of  these  equations  is  quite  similar  to  that  of  [4] — i.e., 
each  estimator  consists  of  an  infinite  band  of  filters,  and  the  filter  for  the 
p moment  is  coupled  only  to  those  for  the  first  and  (p+1)  moments.  There- 
fore, we  are  led  to  the  design  of  suboptlmal  estimators.  The  technique  pro- 
posed here  is  motivated  by  the  highly  successful  use  of  folded  normal  assumed 
density  approximations  in  the  phase  tracking  problem  [4];  filters  designed 
using  this  technique  performed  very  well  as  compared  with  other  suboptlmal 
estimators.  We  will  describe  similar  techniques  for  the  design  of  suboptlmal 
estimators  on  s”  and  SO(n) . ' 


We  first  review  the  notions  of  Brownian  notion  and  Gaussian  densities 
on  Lie  t^roups  and  ho?!Ogeneous  spaces.  Yosida  [28]  proved  that  the  funda^ 
mental  solution  of 

..  -yApCx,.)  = 0 (20) 

where  y > 0 and  A is  the  Laplace-Beltrami  operator  CLaplaclan)  on  a Rieman- 
nian  homogeneous  space  M [5,13,29],  is  the  density  (with  respect  to  the  Rie- 
mannlan  measure)  of  a Brownian  motion  on  M 

According  to  [5],  the  fundamental  solution  of  (20)  is  given  by 

(t~tg)y 

pCx,t;xQ,tQ)  = Z <}>^Cx)<j>^CxQ)e  (21) 

/•  i 

v\?here  and  cp.  are  the  eigenvalues  and  the  corresponding  eigenfunctions  of 
the  Laplacian  (see  the  Appendix).  The  function  p(x,t;  ^q/^q)  the  solution 
to  (20)  v/ith  initial  condition  equal  to  the  singular  distribution  concentrated 
at  X = Xq.  Also,Grenander  [27]  defines  a Gaussian  (normal)  density  to  be  the 
solution  of  (20)  for  some  t. 

• The  folded  normal  density  F(9;T1,y)  used  by  Willsky  as  an  assumed 
density  approximation  for  the  phase  tracking  problem  is  indeed  a normal 
density  on  in  the  sense  of  Grenander  [4];  in  fact,  the  trigonometric 
polynomials  e are  eigenfunctions  of  the  Laplacian  on  S . Motivated 

by  the  success  of  Willsky's  suboptimal  filter,  vre  will  design  subopUmal 
esUmators  for  the  SO(n)  and  bilinear  systems  by  employing  normal  assumed 


^ Yosida  defines  a Brownian  motion  process  to  be  a temporally  and  spatially 
homogeneous  Markov  process  on  M which  satisfies  a continuity  condition  of 
Llndeberg's  type. 
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^&£ 


■1  n 


conditional  densities  of  the  form 


-X.YCt) 

p(x,t)  ==  Z (f>^(x)(|!^Cn(t))e  ^ 


where  n(,t)  and  Y(t)  are  parameters  of  the  density  which  are  to  be  estimated,^ 


^ In  order  to  assure  the  existence  of  a conditional  density,  it  Is  sufficient 
to  assume  that  the  system  is  "controllable  from  the  noise"  [10,13,17J. 


In  this  section  ^he  suboptimal  estimation  technique  discussed  in 

the  previous  section  will  be  used  in  order  to  design  filters  for  the  S*' 

estimation  problem  Cl),C9),  The  optimal  constrained  least-squares  estimator 

is  described  by  (14)  and  (18)- (19).  First,  the  suboptimal  estimator  for 
2 

S will  be  described  in  detail;  then  the  generalization  to  s“  will  be 
discussed, 

2 

In  our  discussion  of  estimation  on  S , we  will  refer  to  a point  on 

2 A 

S in  terms  of  the  Cartesian  coordinates  x = (Xj^,X2,X2)  or  the  polar 

coordinates  (d,4>)  (see  the  Appendix,  in  which  harmonic  analysis  on 

is  summarized) . According  to  the  Appendix,  the  spherical  harmonics 

(defined  in  (A. 11)- (A- 12))  are  the  eigenfunctions  of  the  Laplacian  A £ 

s 

(defined  in  (A. 10)),  and  all  spherical  harmonics  of  degree  i have  the  same 

eigenvalue  -1(Z+1).  Thus  the  assumed  density  approximation  is  a normal 
2 

density  on  S of  the  form  (22) , as  discussed  in  the  previous  section: 

. p(e,cp,t)=  £ £ Y (e,cp)Y*  (23) 

£=0  m=-£ 


where  * denotes  complex  conjugate.  In  other  words. 


“ E’'[Y*^(e(t),Kt)]  = y*J.n(t).A(t))e'^^‘^'^^^^^*'^  . (24) 


In  order  to  truncate  the  optimal  estimator  (18)-(19)  after  the 
x^^^(t|t)  equation  using  the  assumed  density  (.23),  we  must  compute 

(t)x' (t)  ] , or  equivalently,  x^^^^^(t|t),  in  terms  of  x^^^(t|t), 
p = 1,2,...,N.  However,  if  x(t|t)  is  known,  so  are  c ^(t)  and  c^,(t).  A 


11 


simple  computation  1133,  [22]  then  shows  that  ^,m  = -^(N+l) ,,,,  ,N+l) 


can  be  corap u ted  from 


N+: 


= (-1) 


m 


. [(N+l-m)! 

2N+3l2  ^ / 

o 

o 

[(N+  1+m)! 

4nJ 

/ 

wr/'  r. 

/ o 9\  ~] 

i(N+l){N+2) 

rj  B 

r (^10^^^'^  ^ 1 )j 

(25) 

where  P,  are  the  associated  Legendre  functions  (see  (A, 11)- (A- 12) ) , Finally, 
£m 


it  is  shown  in  [13]  that  there  exists  a nonsingular  matrix  T such  that 

fN+l]  , 


Tx‘ 


.m-ii 


(26) 


Ihu=  t)  can  be  computed  from  -(N+1)  < m < N+l)  and 


x^^~l^(t[t)-  The  optimal  estimator  (18)  is  truncated  by  substituting  this 


approximation  for  |t)  into  the  equation  for  x^^^(t  |t) . Notice  that 


the  entire  procedure  for  truncating  the  optimal  estimator  can  equivalently 
be  performed  on  the  infinite  set  of  coupled  equations  for  the  generalized 


Fourier  coefficients  c.  (t),  using  the  approximation  (24) 

£m 


We  note  that  one  can  shov/  that 
Of  (t)  = yi!x(t  |t)  II  ^ 1 


and  this  quantity  can  be  used  as  a measure  of  our  confidence  in  our  estimate. 
If  x(t[t)  satisfies  the  assumed  density  (23)  , 

(27) 


«(t)=  |!x(t|t)l|  = e'Y(^) 


so  Y = 0 (zero  "variance  ')  implies  or  = 1 , and  y = “ (infinite  "variance")  Implies 
a = 0 (see  [4]  for  the  analog). 


!i 


: i. 
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2 

Example  1;  Suppose  that  we  truncate  the  optimal  S estimator  (18) 

■»  [21  , 

after  N =>  1— i,e.,  we  approximate  (tit)  using  the  above  approximation. 

Assume  that  Q(t)  = I and  {A^,  i= 1,2,3)  are  given-  Then  the  resulting 
suboptimal  estimator  is  (for  Q(.t)  = I) 

Clx(t|t)  = [A  + ^ E A^]  x(t  |t)dt 
^ i=l 


+ P(t)H'{t)R“^t)[dz2(t)  - H(t)x(tjt)dt]  (28) 

where  the  "covariance"  matrix  P(t)  is  given  by 

Pjpt)  =ij(t|t)(|  |ii(t|t)||  - 1)  (29) 

for  i j , j 5^  k , and 

PijW  =Xj(t)t)Xj(tlt)(llx(tlt)ll  -1)  (30) 

for  i 3*  j.  It  is  shown  in  [37]  that  the  matrix  P(t)  of  (29)-(30)  is  positive 
semideflnite,  and  thus  can  be  viewed  as  a covariance  matrix.  The  results  .of 
Monte  Carlo  simulations  to  evaluate  the  performance  of  this  estimator  are 
presented  in  Section  VII. 

• The  extension  to  S*'  of  -his  technique  for  constructing  suboptimal  • 
estimators  is  straightforward.  The  procedure  uses  the  spherical  harmonics_ 

on  S*'.  In  polar  coordinates,  a point  on  S^.can  be  described  by  - -■  — 
^®1’®2’ " ‘ ’-n-l’*^^  = (0,i{>),  where  0 < 0^  < ir  and  0 < ({» -<  2Tr.  Also,  the 
spherical  harmonics  are  denoted  by 




n-2 

n (sine  ) 

k=0 


c 


m,  -m,  , 
k k+1 


(cosSk^P 


..VT->4.r 


Li 


where  i s m,  s . . . sm  , ^ 0 and  are  the  Gegenbauer  polynomials  [33] 
1 n-1  j 


(that  is,  the  functions  are  eigenfunctions  of  the  Laplace-Bei  trami 


operator  with  eigenvalue  -il(n+£-l)).  Hence  the  assumed  density  approximation 


on  S iS 


P(6,,,l)=  S (32, 

I,  im; 


That  is,  c„  , .(t)  :4e\y*  , . (e(t) , m(t))]  is  assumed  to  be 


-I  (i+n-l)Y(t) 


(33) 


The  procedure  for  truncating  the  filter (18)  is  identical  to  the  case. 

If  x(t|t)  is  known,  so  are  c . .(t) , and  these  can  be  used  to  compute  Y(t) , 

^ ArTj4- 1 1 

I](t),  and  \{t).  Then  [c..  . } can  be  computed  from  (33),  and  ^"^(tlt) 

can  be  computed  from  ^ estimator  is  truncated 


"[N+1], 


by  substituting  this  approximate  expression  for  x^  ■^(t|t)  into  the  equation 


(18)  for  x’-^^(t|t). 


lA 


: i 


IV.  ESTIMATION  ON  SO(n) 


In  this  section  we  discuss  the  construction  of  suboptimal  estimators 


for  the  SO(n)  estimation  problem  (7)-(8).  We  will  only  consider  the  S0(3) 
problem;  the  results  are  extended  to  SO(n)  in  [13].  T!ie  concepts  of  harmonic 
analysis  on  S0(3)  presented  in  the  Appendix  will  be  used  extensively. 

I 

Consider  the  sequence  {D  , f = 0, 1 , . . . } of  irreducible  unitary 

representations  of  SO(3) , as  defined  in  (A.5)-(A.6).  Theorem  A.  1 implies 

that,  for  fixed  I , the  matrix  elements  -fsm,nif}  are  eigenfunctions 

of  the  Laplacian  (defined  in  (A. 4)) with  the  same  eigenvalue  X ; also, 

all  eigenfunctions  of  the  Laplacian  can  be  written  as  linear  combination  of 

the  fD"^  1.  Hence,  the  assumed  density  which  will  be  used  to  truncate  the 
'■  mn  ^ 

optimal  estimator  (16)-(17)  is  a normal  density  on  SO(3)  of  the  form  (22): 

p(R.t)=  S i e '(34) 


1=0  m,n=-f 


where  R,  Tl(t)  e SO(3)  and  Y(t)  is  a scalar.  That  is, 

cLw  = E^[D;J,^(71(t))*3 

mn  mn 


is  assumed  to  be 


-X,Y(t) 


cf,  (t)  = (TlCt))*  e ^ . (36) 

mn  mn 

The  procedure  for  truncating  the  filter  (16)  is  similar  to  the  case, 

although  we  make  use  of  some  additional  concepts  from  representation  theory. 

* 1 1 
If  X(t)t)  is  known,  so  are  -1  s m.n  s 1},  since  D is  equivalent  to 

the  self-representation  of  SO(3).  Define  the  matrix  C^(t)  v/ith  elements 

c"*  (t) , -X  s:  m,n  s jj;  then 
mn 

A -1  1 1 • 1 * ■•2XiY(t) 

A(t)  A C^t)cNt)  = [D^-n(t))]  [D^C-nCt))]  e • 
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-2X,Y(t) 

= I . e ^ (37) 

since  is  unitary  (here  C is  the  hermitian  transpose  of  C) , Thus  Y(t)  can 
be  computed  from 

Y(t)  = log  [i  tr  A(t)]  . (38) 

Then  the  elements  of  Tl(t)  can  be  computed  from  (36)  and  (38),  since  D^(Tl(t)) 

N+1 

is  similar  to  Ti(t) . Once  Y(t)  and  7l(t)  have  been  computed,  {c^^  ; 

“(N+l)  ^m,n  tsN+l]  are  computed  from  the  formula  (36). 

In  order  to  truncate (16)  after  the  moment  equation,  we  must 

approximate  E^[X^^^(t)h'(t)X(t)];  however,  this  matrix  consists  of  time-varying 
^ “tN+l] 

deterministic  functions  multiplying  elements  of  X (t  jt),  so  we  will  show 
how  to  approximate  this  matrix.  The  symmetrized  Kronecker  p^h  power 
operating  on  the  symmetric  tensors  x^^^such  that  = 1 furnishes 

a representation  of  SO(3)  which  is  reducible  [26].  It  is  shown  in  [13]  that 
there  is  a nonsingular  matrix  T such  that 


D^(X) 


.[p-23 


The  matrix  T is  related  to  the  Clebsch-Gordan  coefficients  [6] , but  T can 

also  be  computed  by  the  method  of  Gantmacher  [34,  p.  160].  It  is  clear 

from  the  decomposition  (39)  that  can  be  computed  from  C^^^(t) 

and  X^^  ^^(t|t).  The  optimal  estimator  (16)  is  truncated  by  substituting 

this  approximation  into  the  equation  for  X {t|t) . 

We  note  here  that,  due  to  the  decomposition  (39),  the  estimation 

equations  and  the  truncation  procedure  could  have  been  expressed  solely 

in  terms  of  the  irreducible  representations  D^(X(t)).  However,  we  have 

chosen  to  work  with  the  X^^^  equations  primarily  for  ease  of  notation.  For 

large  N,  the  equations  would  provide  significant  computational  savings 
Fnl 

over  the  equations,  as  these  are  redundant;  however,  the  practical  , 
implementation  of  this  technique  will  probably  bo  limited  to  small  values  of  N. 


V.  SIMULATION  RESULTS 


As  an  illustration  of  the  techniques  presented  in  the  previous  sections, 
the  first  order  filter  (FOF)  of  Example  1 (Section  III)  was  evaluated  by  means 
of  digital  Monte  Carlo  simulations.  It  was  compared  to  both  the  extended 
Kalman  filter  (EKF)  [7]  and  the  Gustafson-Speyer  linear,  minimum-variance 
quadrature  filter  (LQF)  [24],  Identical  noise  sequences  v/ere  used  to  allow 
direct  comparisons . 

2 

The  system  considered  was  the  S system,  i.e., 

3 

dx(t)  = Fx(t)  dt -f-  E A.  x(t)dw  (t)  (40) 


1/2 

dz(t)  = x(t)  dt  + r ' dv(t) 

3 1 3 2 

where  F = E f .A. +-r-q  E A.  , and  {A,  ,1=1,2, 3}  are  the  skew-symmetric 
i=l  ^ ^ i=l  ^ 


(41) 


matrices 


0 0 0 

0 0-1 

0 10 


0 0 1 

0 0 0 

-10  0 


0-10 
10  0 
0 0 0 


(42) 


Also,  w(t)  has  strength  ql;  and  v has  strength  I.  In  this  experiment. 


the  nominal  angular  velocities  {f^,i=l,2,3l  w’ere  chosen  to  be  100.0,  and 
q and ■ r were  varied . 

rw  • ^ ‘ tk  - 

For  all  three  filters,  the  normalized  estimate  x(t)  = x(t|t)/||x(t|t)||  was 
used.  The  filters  have  an  identical  structure  for  the  approximate  x equation: 

didjt)  = Fx(t|t)  + ^P(t)  [dz(t)  - x(t|t)dt]  . (42) 

However,  for  the  FOF,  P(t)  is  given  by  the  highly  nonlinear  memoryless 
equations  (29)-(30).  In  the  EKF,  P(t)  satisfies  the  Riccati  equation 


^ P(t)  = FP(t)  + P(t)F*  + qGCktIt))e’  (x(t|t)  - P(t)P'  (t)  (43) 

where 

G(x)  = [AjX,  A2X,  Agx]  . (44) 

$ 

Since  the  Riccati  equation  (43)  is  a function  of  x(tlt),  the  P(t)  calculation 
in  the  EKF  requires  extensive  on-line  computation,  which  represents  a con- 
siderable burden.  17 


In  the  LQF,  P(t)  is  given  as  the  solution  of  the  coupled  dynamic  equations 


^ P(t)  = FP(t)  + P(t)F’  + A(X(t),t)  - - P(t)P'(t) 
dC  r 


^ X(t)  = FX(t)  + X(t)F'  + A(X(t).t) 
where  A(X(t),t)  is  a diagonal  matrix  with  1th  component 


MX(t).t)  - , I (A„)  (47) 

k,£,ra=l 

Notice  that  these  equations  for  P(t)  and  X(t)  can  be  calculated  off-line, 
but  the  LQF  thus  has  a considerable  storage  requirement.  Because  P(t)  in 
the  FOF  is  only  given  by  a memoryless  nonlinearity,  this  filter  requires 
considerably  less  storage  than  the  LQF  and  less  on-line  computation  than  the 
EKF. 

Our  approach  to  the  statistical  analysis  of  the  Monte  Carlo  simulations 
closely  parallels  that  of  Bucy  and  his  associates  [l],  [3],  [35],  The  steady- 
state  mean-squared  error 

3 

' P-  = E[llx(t)  - mf]  = S E[(x  (t)  - x(t))h  (48) 

^ 1=1  ^ ^ 

where  x^(t)  denotes  the  estimate  of  the  ith  component  of  the  state  x^Ct),  was 

n 

used  as  the  performance  criterion.  If  {x  } and  {x  },  i = 1, . . . ,N,  are  se- 
quences of  independent  realizations  of  x(t)  and  5l^(t) , respectively,  then  the 


statistic 


n ~nii2 


is  an  approximation  to  sufficiently  large  N.  In  fact,  by  the  Central 

Limit  Theorem  [36,  p.  278],  ^2  is  asymptotically  normal  with 

E[{l2]  = (50) 

var  [^2]  = ^ + 2(m-4)j2  + + ^(^4^23  “ ^^*2^  ( (51) 

i 

j.  $ 

where  = E[(x^(t)  - x^(t))  ] and 


<. 


mi  ^ 


^ E[(x^(t)  ~ x^{t))^(x^(t)  - Xj(t))^]. 

Thus,  for  large  N,  with  probability  0.9974,  the  3ct  confidence  inter\.'al  is 
given  by 

I ^2  ~ ^2  I ^ ^ v/var  (C^)  » (52) 


or  equivalently. 


Pr 


^2  ^ ^2 
1 + 3^/a  ^2  ^ 1 - 3Va 


= 0.9974 


where 


a = 


var  (n-2) 

(^2)' 


(53) 


(54) 


In  the  Monte  Carlo  simulations , a was  estimated  from  the  samples  (using 
sample  means  as  in  (49)),  and  approximate  confidence  intervals  were  thus 
computed . 

In  the  experiment,  15  sample  paths,  each  of  which  contained  1000 
steps  of  length  ,001  seconds,  were  run  in  each  simulation.  The  first  200 
samples  in  each  sample  path  were  discarded  to  allow  the  transients  to  decay, 
so  the  remaining  800  samples  represented  steady-state.  If  all  the  steady- 
state  errors  were  averaged  as  in  (49),  this  would  lead  to  12000  samples  of 
the  steady-state  error.  However,  as  noted  in  [4],  [24]  and  [35],  adjacent 
errors  in  each  path  are  correlated,  so  the  effective  Monte  Carlo  length  is 
somewhere  In  the  range  between  N=1200  and  N=12000.  The  three  standard 
deviation  confidence  intervals  were  calculated  for  both  values  of  N. 

The  results  of  the  simulations  are  presented  in  Table  I.  The  3a 
confidence  intervals  Ij^  (for  N = 12000)  and  I2  (for  N = 1200)  are  shown. 

The  results  of  this  approximate  statistical  analysis  of  the  Monte  Carlo 

simulations  indicate  that,  for  this  simple  example,  the  FOF 
performs  comparably  to  the  l.QF,  and  better  than  the  EKF . , 

The  FOF  seems  to  perform  better  in  comparison  to 
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1 


the  other  filters  as  q increases,  due  to  the  increasing  dominance  of  the 
bilinear  noise  term  in  the  system  equation  (AO).  These  results  are  significant, 
due  to  the  fact  that  the  FOF  designed  here  requires  considerably  less  storage 
and  computation  than  the  other  filters  (no  additional  differential  equations  or 
storage  for  P(t)  are  required). 


M.S.  Error  Ig 

U2  (N=12000)  (N=1200) 


q 

= 0.01 

FOF 

9.65 

[9.41,9.91] 

[8.92,10.52] 

= 0.01 

EKF 

10.21 

[9.96,10.48] 

[9.45,11.11] 

r 

LQF 

9.49 

[9.26,9.75] 

[8.78,10.34] 

q 

= 0.01. 

FOF 

11.69 

[11.39,12.01] 

[10.7^,12.77] 

= 1.00 

EKF 

11.75 

[11.44,12.07] 

[10.84,12.82], 

r 

LQF 

11.69 

[11.39,12.01] 

[10.79,12.77] 

q 

= 1.00 

FOF 

170.38 

[165.26,175.82] 

[155.19,188.87] 

= 0.01 

EKF 

193.75 

[188.74,199.04] 

[178.74,211.511 

r 

LQF 

170.65 

[165.58,176.03] 

[155.60,188.92] 

Table  I.  Monte  Carlo  M.S.  Estimation  Error(s) 

(xlO"^) 


$ 
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VIII.  CONCLUSIONS 

The  state  estimation  problem  for  bilinear  stochastic  systems  evolving 
on  compact  Lie  groups  unci  homogeneous  space  has  been  consiciered.  The 
techniques  of  harmonic  analysis  on  compact  Lie  groups  have  been  applied  to 
the  design  of  suboptimal  estimators  for  such  systems.  Monte  Carlo  simula- 
tions of  a simple  example  indicate  that  a computationally  simple  filter  designed 
by  these  methods  performs  favorably  as  compared  to  two  other  filters. 


■j 


i 
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APPEMDIX  - HARMONIC  ANALYSIS  ON  CO>iPACT  LIE  GROUPS 


In  this  appendix  we  suuiiiiarize  socie  facts  fioji  the  theory  of  harmonic 
analysis  on  Lie  groups.  For  details  of  group  representations  see  references 
15,6,13,29,301. 

Definition  A.l:  A flnlte-dlmensional  natrix  representation  of  a 
compact  Lie  group  G Is  a continuous  homomorphism  D which  maps  G into  the 
group  of  nonsingular  linear  transformations  on  a finite  dimensional  vector 
space  V.  The  representative  ring  of  G is  the  ring  generated  over  the  field 
of  complex  numbers  by  the  set  of  all  continuous  functions  on  G which 
are  matrix  elements  of  some  irreducible  representation  D. 


The  Peter-Weyl  theorem  is  the  major  result  in  harmonic  analysis  on 
compact  Lie  groups  [5,6,29,30].  It  gives  the  direct  sum  decomposition  of 
L2(G),  the  space  of  square  integrable  functions  with  respect  to  the  Haar 


measure  dg: 


L„(G)  = © , ^ 

^2'  aeA  a 


(A.l) 


where  A is  the  set  of  equivalence  classes  of  finite  dimensional  irreducible 
representations  of  G,  and  denotes  the  vector  space  spanned  by  the  (n^)*" 
functions  fD®. ; ij=l , . . . ,n  3. 

The  Laplace-Beltrami  operator  (Laplacian)  on  a Lie  group  or  Riemannian 
homogeneous  space  is  discussed  in  [S],  [13]  and  [29].  Some  examples  will  be 
presented  here,  but  one  very  important  result  for  our  purposes  is  the  following. 
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Theorem  A.  1 [5,  p.  40],  [29,  p.  257];  Let  G be  a compact  Lie  group, 

and  let  V be  defined  as  in  equation  (A.l).  Then  each  function  cp  e V is  an 
a a 

eigenfunction  of  the  Laplacian  A/  ^tid  all  o e V have  the  same  eigenvalue 

f/ 

X . Conversely,  each  eigenfunction  9 of  the  Laplacian  is  an  element  of  the 
a 

representative  ring. 

Hence,  harmonic  analysis  on  a compact  Lie  group  can  be  performed 
either  in  terms  of  the  representative  ring  or  the  eigenfunctions  of  the 
Laplacian,  since  these  two  sets  of  functions  are  the  same.  In  this  paper  we 
are  primarily  concerned  with  the  application  of  these  results  to  the  special 
orthogonal  group  SO(n)  and  the  n-sphere 
The  Lie  group  SO(n)  is  defined  by 


SO(n)  = {X  eIR'^^"lX'X=  I,  detX  = + l}  . 


(A.2) 


The  theory  of  representations  of  SO(n)  is  discussed  in  [13],  [26];  for  this 
paper  we  need  only  consider  SO(3) . Any  matrix  R in  SO(3)  can  be  described 
in  local  coordinates  in  terms  of  the  Euler  angles  9,  6,  ip,  domain 
0 < <{>  < 2ir,  O<0<2Tr,O<ip<2Tr  [6].  An  element  of  S0(3)  will  thus 

be  denoted  by  R(<J>,0,'P)  or  just  ((j),0,ip). 

In  the  Euler  angle  coordinates,  the  (unnormalized)  Haar  measure  is 


. dij,(cp,  e,  Y)  = sinBdcpdedY 

and  the  corresponding  Laplace-Beltrami  operator  is  given  by  [3  1] 


(A.  3) 


2 2 2 

^so(3)  " 7i^  ae  ~ ^ 

Talman  [6]  computes  a sequence  D '(cp,  0,  Y) , / = 0, 1 , . . . , of  unitary 
irreducible  representations  of  SO(3);  its  matrix  elements  are  given  by 


d‘  fa,e,n  = r-"e-“"‘fd‘  {e)e- 

mn  mn 


(A.  5) 


v/here 
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for  -X  ^ m,n  s X.  Here  t is  summed  over  all  nonnegative  integers  such  that 
the  arguments  of  the  factorial  functions  in  (A. 6)  are  nonnegative;  i.e.. 


m-nstsx+m,  Ostsx-n  . 

In  fact,  these  are  (up  to  equivalence)  all  of  the  irreducible  representations  of 
SO(3).  The  Peter-Weyl  Theorem  yields  the  decomposition 

L^^OCa))  = © V (A, 7) 

X ^ 

where  is  the  vector  space  spanned  by  the  (2x+l)^  functions 
f^mn'  

The  n-sphere  S*'  = = 1 } is  dif feomorphic  to  the  homoge- 

neous space  SO(n) /SO(n-l) . Harmonic  analysis  on  is  studied  in  terms  of 
the  spherical  harmonics  [6,30,33].  The  space  of  spherical  harmonics  of 
degree  X can  be  characterized  as  the  eigenspace  of  the  SO (n+1) -invariant 

A _ on  S*'  with  eigenvalue  -X(n-l+£);  other  equivalent  characterizations  are 
S 

given  in  [13]. 

2 

In  particular,  we  consider  the  2-sphere  S . Any  point  (Xj^,X2,x^) 

2 

on  S can  be  expressed  in  the  polar  coordinates  (0,<(>)»  where  0 < 0 < w, 

0 < it  < 2ir,  by  defining 

Xj  = cose;  = sine  coscp;  x^  = sine  sincp  . (A. 8) 

• The  Riemannian  measure  invariant  under  the  action  of  SO(3)  is 

d^l{e,cp)  = sin9  d edcp  . (A. 9) 

I 

The  corresponding  invariant  Laplace-Beltrami  operator  is 


2A 


Y (a  \ - riLlIBlL  (2l±. 

- <-*>  [(£+m)!  4n 


J In 


(cos  9) 


(A. 11)  ' 


(A. 12) 


for  jj=0, 1 , , . . and  m = 0 .1 .....  I,  where  P (cos  9)  are  the  associated 

i m 


Legendre  functions  and  * denotes  complex  conjugate.  Notice  that  Y.  Is  an 

XrQ 


eigenfunction  of  A „ with  eigenvalue  . 

S 


REFERENCES 

1.  R.S,  Bucy,  C.  Hecht,  andK.D.  Senno,  ’'New  Melhocis  for  Nonlinear 
Filtering/'  Revue  Frcincaise  ci'Aulomatique , Informalique  et  cie 
Recherche  Operationnelle , February  1973,  pp.  3-54. 

2.  R.S.  Bucy  and  H.  Youssef,  "Fourier  Realization  of  the  Optimal  Phase 
Demodulator,"  Fourth  Symposium  on  Nonlinear  Estimation  and  Its 
Applications,  San  Diego,  CA.  , 1973. 

3.  J.I.  Mallinckrodt,  R.S.  Bucy,  and  S.Y.  Chang,  Final  Project  Report 
for  a Design  Study  for  an  Optimal  Non-Linear  Receiver/Demodulator, 

NASA  Contract  NAS  5-10789,  Goddard  Space  Flight  Center , Maryland, 
1970. 

4.  A.S.  Wlllsky,  "Fourier  Series  and  Estimation  on  the  Circle  with 
Applications  to  Synchronous  Communication,  Parts  I,  II,"  IEEE  Trans . 
on  Inf.  Theory,  Vol,  IT-20,  No.  5,  September  1974,  pp.  577-590. 

5.  E.M.  Stein.  Topics  in  Harmonic  Analysis.  Princeton:  Princeton  Univ. 
Press,  1970. 

6.  J.D.  Talman,  Special  Functions — A Group  Theoretic  Approach.  New 
York:  Benjamin,  1968. 

7.  A.H.  Jazwinski,  Stochastic  Processes  and  Filtering  Theory.  New  York: 
Academic  Press , 1970. 

8.  R.W.  Brocket!,  "System  Theory  on  Group  Manifolds  and  Coset  Spaces," 
SIAM  I.  Control , Vol.  10,  No.  2,  May  1972,  pp.  265-284. 

9.  R.W.  Brocket!,  "Lie  Algebras  and  Lie  Groups  in  Control  Theory," 

in  Geometric  Methods  in  System  Theory,  D.Q.  Mayne  and  R.W.  Brocket!, 
eds.,  Reidel  Publ.  Co.,  The  Netherlands,  1973. 

10.  R.W.  Brocket!,  "Lie  Theory  and  Control  Systems  on  Spheres,"  SIAM 
■ T.  Appl.  Math. . Vol.  25,  No.  2,  Sept.,  1973,  pp.  213-225. 

11.  A.S.  Willsky,  S.I.  Marcus,  and  D.N.  Martin,  "On  the  Stochastic 
Stability  of  Linear  Systems  Containing  Colored  Multiplicative  Noise," 
IEEE  Trans,  on  Aut.  Control.  Vol.  AC-20,  Oct.  1975,  pp.  711-713. 


12.  A.S.  Willsky  and  S.I.  Marcus,  "Analysis  of  Bilinear  Noise  Models 

in  Circuits  and  Devices,"  J.  Franklin  Tnst. , vol.  301,  No.  1,  January 
1976,  pp.  103-122. 


13.  S.I.  Marcus,  Estriraatlon  a»d  Analysis  of  Nonlinear  Stochastic  Systeros, 
Ph.D.  Thesis,  Dept,  of  Electrical  Enginecrlns,  M.I.T.,  Catubridge, 

Mass.,  June  1975;  also  M.I.T.  Electronic  Sys'^ems  Laboratory  Report 
No.  ESL-R-601. 

14.  S.I.  Marcus  and  A.S.  Willsky,  "Algebraic  Structure  and  Finite  Dimen- 
sional Nonlinear  Estimation,"  SIAM  Journal,  on  Mathematical  Analysis, 
to  appear. 

15.  J.T.  Lo,  "Signal  Detection  on  Lie  Groups,"  in  Geometric  Methods  in 
System  Theory,  D.Q.  Mayne  and  R.W.  Brockett,  eds.,  Reldcl  Pub.  Co., 

The  Netherlands,  1973. 

16.  J.M.C.  Clark,  "An  Introduction  to  Stochastic  Differential  Equations 
on  Manifolds,"  in  Geometric  Methods  in  System  Theory,  D.Q.  Mayne  and 
R.W.  Brockett,  eds.,  Reidel  Pub.  Co.,  The  Netherlands,  1973. 

17.  D.L.  Elliott,  "Diffusions  on  Manifolds  Arising  from  Controllable  Sys- 
tems," in  Geometric  Methods  in  System  Theory,  D.Q.  Mayne  and  R.W. 
Brockett,  eds.,  Reidel  Pub.  Co.,  The  Netherlands,  1973. 

18.  J.T.H.  Lo,  "Exponential  Fourier  Densities  and  Optimal  Estimation  and 
Detection  on  the  Circle,"  IEEE  Trans.  Inform.  Theory,  vol.  IT-23, 

January  1977,  pp.  110-116. 

2 

19.  J.T.H.  Lo  and  L.R.  Eshleman,  "Exponential  Fourier  Densities  on  S and 
Optimal  Estimation  and  Detection  for  Directional  Processes,"  IEEE 
Trans.  Inform.  Theory,  Vol.  IT-23,  May  1977,  pp.  321-336. 

20.  J.T.H.  Lo  and  L.R.  Eshleman,  "Exponential  Fourier  Densities  on  S0(3) 
and  Optimal  Estimation  and  Detection  for  Rotation  Processes,"  in 
Proc.  1976  IEEE  Conf.  on  Decision  and  Control,  pp.  134-142. 

21.  A.S.  Willsky  and  S.I.  Marcus,  "Estimation  for  Bilinear  Stochastic 
Systems,"  in  Variable  Structure  Systems  with  Applications  to  Economics 
and  Biology,  A.  Ruberti  and  R.R.  Kohler,  eds..  New  York,  Springer- 
Verlag,  1975,  pp.  116-137. 

22.  S.I.  Marcus  and  A.S.  Willsky,  "The  Use  of  Spherical  Harmonics  in  Sub- 
optimal  Estimator  Design,"  in  Proc.  1976  IEEE  Conf.  on  Decision  and 
Control,  pp.  199-203. 

23.  I.Y.  Bar-Itzhack,  "Iterative  Optimal  Orthogonallzation  of  the  Strap- 
down  Matrix,"  IEEE  Trans.  Aero,  and  Elec.  Systems,  Vol.  AES-11,  January 

1975,  pp.  30-37. 

24.  D.E.  Gustafson  and  J.L.  Speyer,  "Linear  Minimum  Variance  Filters  Applied 
to  Carrier  Tracking,"  IEEE  Trans.  Aut.  Control,  Vol.  AC-21,  February 

1976,  pp.  65-73. 

$ 

25.  G.L.  Blankenship,  "Lie  Theory  and  the  Motnent  Stability  Problem  in  Sto- 
chastic Differential  Equations,"  Triennial  World  Congress  of  IFAC, 
Cambridge,  MA.,  August  24-30,  1975. 

27 


26.  F.D.  Murnaghan,  Theory  of  Group  Representations.  Baltimore: 

Johns  Hopkins  Press,  1938. 

27.  U.  Grenander,  Probabilities  on  Algebraic  Structures.  New  York: 

John  V/iley,  1963. 

28.  K.  Yosida,  Functional  Analysis.  New  York:  Springer-Verlag,  1974, 
pp.  398-403.' 

29.  F.  Warner,  Foundations  of  Differential  Manifolds  and  Lie  Groups. 

Glenview,  111.:  Scott,  Foresman  and  Co. , 1971. 

30.  C.F.  Dunkl  and  D.E.  Ramirez,  Topics  in  Harmonic  Analysis.  New 
York:  Appleton-Century-Crofts , 1971. 

31.  R.  Burridge,  G.C.  Papanicolaou,  and  D.  McLaughlin,  "The  Solution  by 
Descent  of  Certain  Partial  Differential  Equations  on  the  Sphere  and 
Lobachevsky  Plane,"  Comm.  PureAppl.  Math.,  Vol.  XXVI,  1973,  pp. 
105-129. 

32.  R.W.  Brockett,  Finite  Dimensional  Linear  Systems.  New  York:  John 
Wiley,  1970. 

33.  A.  Erdelyi,  W.  Magnus,  F.  Oberhettinger,  and  F.  Tricomi,  Higher 
Transcendental  Functions,  Vol.  II.  New  York:  McGraw-Hill,  1953. 

34.  F.R.  Gantmacher,  The  Theory  of  Matrices,  Vol.  I.  New  York:  Chelsea, 
1960. 

35.  R.S.  Bucy,  C.  Hecht,  and  K.D.  Senne,  An  Engineers*  Guide  to  Building 
Nonlinear  Filters,  F.J.  Seiler  Research  Lab.  Report  No.  SRL-TR-72-0004, 
May  1972. 

36.  H.  Bauer,  Probability  Theory  and  Elements  of  Measure  Theory.  New 
York:  Holt,  Rinehart,  and  Winston,  1972. 

37.  K.  Hsu,  The  Evaluation  of  Several  Suboptlmal  Filters  for  Bilinear 
Systems  Evolving  on  Spheres,  M.S.  Thesis,  Dept,  of  Electrical  Engineering, 
University  of  Texas  at  Austin,  Austin,  Texas,  Kay  1977. 


1 

i i 
f 


28 


